Variation in vocal behavior under varying stress levels
Budgie INBRE stress and contact calls
Marcelo Araya-Salas, PhD & Timothy Wright, PhD
24-02-2022
Load packages
## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2", "kableExtra", "knitr", "formatR", "randomForest", "Rtsne", "MASS", "adehabitatHR", "sp", "GGally", "spatstat", "raster", "vegan", "maRce10/PhenotypeSpace")
aa <- lapply(x, function(y) {
# get pakage name
pkg <- strsplit(y, "/")[[1]]
pkg <- pkg[length(pkg)]
# check if installed, if not then install
if (!pkg %in% installed.packages()[,"Package"]) {
if (grepl("/", y)) remotes::install_github(y, force = TRUE) else
install.packages(y)
}
# load package
a <- try(require(pkg, character.only = T), silent = T)
if (!a) remove.packages(pkg)
})Functions and global parameters
opts_knit$set(root.dir = "..")
# set evaluation false
opts_chunk$set(fig.width = 6, fig.height = 3, warning = FALSE, message = FALSE, tidy = TRUE)
read_excel_df <- function(...) data.frame(read_excel(...))
# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")
# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd", "freq.median",
"freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
"mindom", "maxdom", "dfrange", "modindx", "meanpeakf")]
sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])
sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year, sep = "-")
sels$date <- as.Date(sels$date, format = "%d-%b-%Y")Counts per individual
df <- as.data.frame(table(sels$ID))
names(df) <- c("ID", "Sample_size")
df <- df[order(df$Sample_size, decreasing = FALSE), ]
kb <- kable(df, row.names = FALSE)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)| ID | Sample_size |
|---|---|
| 132YGMM | 6 |
| 125YGMM | 12 |
| 160YGHM | 16 |
| 310YGHM | 21 |
| 393YGHM | 25 |
| 303YGHM | 37 |
| 398WBHM | 41 |
| 365YGHM | 44 |
| 399YGLM | 46 |
| 300YGHM | 47 |
| 270YGHM | 64 |
| 407YGHM | 97 |
| 386WBHM | 100 |
| 377WWLM | 108 |
| 367WWNM | 119 |
| 371YYLM | 149 |
| 397YGHM | 175 |
| 378YYLM | 195 |
| 404WBHM | 196 |
| 258YGHM | 223 |
| 389WWLM | 230 |
| 262YGHM | 306 |
| 400YGHM | 360 |
| 388YGLM | 377 |
| 327YYLM | 404 |
| 396YBHM | 455 |
| 403WBHM | 566 |
| 356WBHM | 761 |
| 361YGHM | 767 |
| 402YGHM | 849 |
| 395WBHM | 854 |
| 360YGHM | 900 |
| 390WBHM | 935 |
| 405YBHM | 975 |
| 385YBHM | 1256 |
| 394YBHM | 1693 |
Add metadata
metadat <- read_excel_df("./data/raw/INBREStress_MasterDataSheet_14Nov19.xlsx")
# head(metadat)
sels$ID[sels$ID == "125YGMM"] <- "125YGHM"
sels$ID[sels$ID == "394YBHM"] <- "394WBHM"
# setdiff(sels$ID, metadat$Bird.ID) setdiff(metadat$Bird.ID, sels$ID)
sels$treatment <- sapply(1:nrow(sels), function(x) {
metadat$Treatment[metadat$Bird.ID == sels$ID[x]][1]
})
sels$treatment.room <- sapply(1:nrow(sels), function(x) {
metadat$Treatment.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$round <- sapply(1:nrow(sels), function(x) {
metadat$Round[metadat$Bird.ID == sels$ID[x]][1]
})
sels$source.room <- sapply(1:nrow(sels), function(x) {
metadat$Source.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$record.group <- sapply(1:nrow(sels), function(x) {
metadat$Record.Group[metadat$Bird.ID == sels$ID[x]][1]
})
# add week
out <- lapply(unique(sels$round), function(x) {
Y <- sels[sels$round == x, ]
min_date <- min(Y$date)
week_limits <- min_date + seq(0, 100, by = 7)
Y$week <- NA
for (i in 2:length(week_limits)) Y$week[Y$date >= week_limits[i - 1] & Y$date <
week_limits[i]] <- i - 1
return(Y)
})
sels <- do.call(rbind, out)
sels$cort.baseline <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$cort.stress <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$weight <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$breath.count <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
# remove week 5
sels <- sels[sels$week != 5, ]Exploratory graph
Physiological parameters
breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count", "D14.Breath.Count",
"D21.Breath.Count", "D28.Breath.Count")])
weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.", "D14.Bird.Weight..g.",
"D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])
cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress", "D14.CORT.Stress",
"D21.CORT.Stress", "D28.CORT.Stress")])
cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline", "D14.CORT.Baseline",
"D21.CORT.Baseline", "D28.CORT.Baseline")])
stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round", "Treatment.Room")],
week = breath.count$ind, breath.count = breath.count$values, weight = weight$values,
cort.stress = cort.stress$values, cort.baseline = cort.baseline$values)
# head(stress)
stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."), "[[", 1),
levels = c("D3", "D7", "D14", "D21", "D28"))
stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)
stress$treatment <- factor(stress$Treatment, levels = c("Control", "Medium Stress",
"High Stress"))
# remove 5th week
stress <- stress[stress$week != "D28", ]
stress_l <- lapply(stress$Bird.ID, function(x) {
X <- stress[stress$Bird.ID == x, ]
X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
X$weight <- X$weight - X$weight[X$week == "D3"]
X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week == "D3"]
return(X)
})
stress <- do.call(rbind, stress_l)
ggplot(data = stress, aes(x = day, y = weight, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Weight (g)", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))
ggplot(data = stress, aes(x = day, y = breath.count, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Breath count", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))
ggplot(data = stress, aes(x = day, y = cort.baseline, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Baseline CORT", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))
ggplot(data = stress, aes(x = day, y = cort.stress, group = Bird.ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Stress CORT", x = "Day sample was taken") +
theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))
Behavioral parameters
Call counts per treatment
call_count <- aggregate(selec ~ week + treatment + round, data = sels, length)
call_count$Week <- factor(call_count$week, levels = (unique(call_count$week))[5:1])
call_count$round <- paste("Round", call_count$round)
call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(data = call_count, aes(x = treatment, y = selec, fill = Week)) + geom_bar(stat = "identity") +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + coord_flip() + facet_wrap(~round,
ncol = 2) + labs(y = "Number of calls", x = "Treatment") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.15))
Call counts per bird and treatment
# Call counts per treatment by individual (in color)
call_count <- aggregate(selec ~ week + treatment + round + ID, data = sels, length)
call_count$Week <- factor(call_count$week, levels = (1:5))
call_count$round <- paste("Round", call_count$round)
call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(data = call_count, aes(x = week, y = selec, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Stress CORT", x = "Day sample was taken") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))Acoustic parameters
call_week_params <- aggregate(cbind(duration, meanfreq, sd, freq.median, freq.IQR,
time.IQR, skew, kurt, sp.ent, time.ent, entropy, meandom, mindom, maxdom, dfrange,
modindx, meanpeakf) ~ week + treatment + round, data = sels, mean)
call_week_params_sd <- aggregate(cbind(duration, meanfreq, sd, freq.median, freq.IQR,
time.IQR, skew, kurt, sp.ent, time.ent, entropy, meandom, mindom, maxdom, dfrange,
modindx, meanpeakf) ~ week + treatment + round, data = sels, sd)
names(call_week_params_sd) <- names(call_week_params) <- c("week", "treatment", "round",
"Duration", "Mean_frequency_(kHz)", "SD_frequency_(kHz)", "Median_frequency_(kHz)",
"Interquantile_frequency_range_(kHz)", "Interquantile_time_range_(s)", "Skewness",
"Kurtosis", "Spectral_entropy", "Time_entropy", "Overall entropy", "Mean_dominant_frequency_(kHz)",
"Minimum_dominant_frequency_(kHz)", "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)",
"Modulation_index", "Mean_peak_frequency_(kHz)")
call_week_params_sd <- call_week_params_sd[, c("Duration", "Mean_frequency_(kHz)",
"SD_frequency_(kHz)", "Median_frequency_(kHz)", "Interquantile_frequency_range_(kHz)",
"Interquantile_time_range_(s)", "Skewness", "Kurtosis", "Spectral_entropy", "Time_entropy",
"Overall entropy", "Mean_dominant_frequency_(kHz)", "Minimum_dominant_frequency_(kHz)",
"Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)", "Modulation_index",
"Mean_peak_frequency_(kHz)")]
names(call_week_params_sd) <- paste(names(call_week_params_sd), "sd", sep = ".")
call_week_params <- cbind(call_week_params, call_week_params_sd)
call_week_params$round <- paste("Round", call_week_params$round)
call_week_params$treatment <- factor(call_week_params$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
pd <- position_dodge(0.3)
ggs <- lapply(c("Duration", "Mean_frequency_(kHz)", "SD_frequency_(kHz)", "Median_frequency_(kHz)",
"Interquantile_frequency_range_(kHz)", "Interquantile_time_range_(s)", "Skewness",
"Kurtosis", "Spectral_entropy", "Time_entropy", "Overall entropy", "Mean_dominant_frequency_(kHz)",
"Minimum_dominant_frequency_(kHz)", "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)",
"Modulation_index", "Mean_peak_frequency_(kHz)"), function(i) {
call_week_params$var <- call_week_params[, i]
call_week_params$var.min <- call_week_params$var - call_week_params[, paste(i,
"sd", sep = ".")]
call_week_params$var.max <- call_week_params$var + call_week_params[, paste(i,
"sd", sep = ".")]
ggplot(call_week_params, aes(x = week, y = var, col = treatment)) + geom_point(size = 4,
position = pd) + geom_errorbar(aes(ymin = var.min, ymax = var.max), width = 0.3,
position = pd, size = 1.7) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = i, x = "Week") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.25))
})
ggs## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
##
## [[12]]
##
## [[13]]
##
## [[14]]
##
## [[15]]
##
## [[16]]
##
## [[17]]
Acoustic space projection
scale_acous_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median",
"freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
"mindom", "maxdom", "dfrange", "modindx", "meanpeakf")])
urfmod <- randomForest(x = scale_acous_param, importance = TRUE, proximity = TRUE,
ntree = 10000)
saveRDS(urfmod, "./data/processed/unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")
mds_urf_prox <- cmdscale(1 - urfmod$proximity, k = 2)
saveRDS(mds_urf_prox, "./data/processed/MDS_unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")Random forest
urfmod <- readRDS("./data/processed/unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")
mds_urf_prox <- readRDS("./data/processed/MDS_unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")
sels$MDS1 <- mds_urf_prox[, 1]
sels$MDS2 <- mds_urf_prox[, 2]
# print(1- sum(urfmod$votes[,1] > urfmod$votes[,2])/ nrow(urfmod$votes))Proportion of real data classified as synthetic data by the unsupervised random forest: 0.00506
ggplot(sels, aes(x = MDS1, y = MDS2, col = as.factor(round))) + geom_point() + labs(color = "Round") +
scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 30) + theme(legend.position = c(0.62,
0.3)) + guides(colour = guide_legend(override.aes = list(size = 10)))
vimp <- as.data.frame(vimp)
vimp$params <- rownames(vimp)
vimp <- vimp[order(vimp$MeanDecreaseAccuracy), ]
vimp$params <- as.factor(vimp$params)
vimp$params <- factor(vimp$params, levels = vimp$params[!duplicated(vimp$params)])
vimp <- vimp[(nrow(vimp) - 15):nrow(vimp), ]
ggplot(vimp, aes(x = MeanDecreaseAccuracy, y = params)) + geom_segment(aes(yend = params),
xend = 0, color = "grey50", size = 1) + geom_point(size = 10, col = viridis(10,
alpha = 0.6)[2]) + labs(x = "Mean decrease accuracy (RF)", y = "Predictors") +
theme(legend.key.size = unit(2, "lines")) + scale_color_discrete(name = "Parameter set") +
scale_fill_discrete(guide = "none") + theme_classic(base_size = 25)
PCA
pca <- prcomp(x = sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
"time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
"maxdom", "dfrange", "modindx", "meanpeakf")], scale. = TRUE)
Y <- as.data.frame(pca$x[, 1:2])
names(Y) <- c("PC1", "PC2")
sels <- data.frame(sels, Y)
ggplot(sels, aes(x = PC1, y = PC2, col = as.factor(round))) + geom_point() + labs(color = "Round") +
scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) + theme(legend.position = c(0.9,
0.8)) + guides(colour = guide_legend(override.aes = list(size = 10)))
t-SNE
scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
"time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
"maxdom", "dfrange", "modindx", "meanpeakf")])
tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE, max_iter = 5000)
saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")
sels <- data.frame(sels, Y)
ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(round))) + geom_point() +
labs(color = "Round") + scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) +
theme(legend.position = c(0.955, 0.25)) + guides(colour = guide_legend(override.aes = list(size = 10)))
Acoustic space by individual (each point) for 3 the dimension reduction methods
- Only individuals with at least 20 calls
- Using rarefaction with n = the minimum sample size across all individuals
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab[tab > 20]), ]
pca_areas <- rarefact_space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID",
parallel = 4, pb = TRUE, type = "mcp")
rf_mds_areas <- rarefact_space_size(X = Y, dimensions = c("MDS1", "MDS2"), group = "ID",
parallel = 4, pb = TRUE, type = "mcp")
tsne_areas <- rarefact_space_size(X = Y, dimensions = c("TSNE1", "TSNE2"), group = "ID",
parallel = 4, pb = TRUE, type = "mcp")
areas <- data.frame(pca_areas[, -ncol(pca_areas)], pca.area = pca_areas$mean.size,
rf.mds.area = rf_mds_areas$mean.size, tsne.area = tsne_areas$mean.size)
saveRDS(areas, "./data/processed/acoustic_space_area_by_individual_for_each_dimension_reduction_method.RDS")areas <- readRDS("./data/processed/acoustic_space_area_by_individual_for_each_dimension_reduction_method.RDS")
ggpairs(areas, columns = which(names(areas) %in% c("pca.area", "rf.mds.area", "tsne.area")),
columnLabels = c("PCA", "Rand.for. MDS", "t-SNE")) + theme_minimal(base_size = 30)
Individual acoustic space across weeks
- PCA used for dimensionality reduction
- Only weeks with at least 6 calls were included
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID.week)
# sort(tab)
Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]
# run with rarefaction
areas_by_week_raref <- rarefact_space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
parallel = 4, pb = TRUE, type = "mcp")
# run without rarefaction
areas_by_week <- space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
cl = 4, pb = TRUE, type = "mcp")
areas_by_week$raref.area <- areas_by_week_raref$mean.size
# add metadata
areas_by_week$ID <- sapply(strsplit(areas_by_week$group, "-"), "[[", 1)
areas_by_week$week <- factor(sapply(strsplit(areas_by_week$group, "-"), "[[", 2),
levels = 1:5)
areas_by_week$treatment <- sapply(1:nrow(areas_by_week), function(x) metadat$Treatment[metadat$Bird.ID ==
areas_by_week$ID[x]][1])
areas_by_week$round <- sapply(1:nrow(areas_by_week), function(x) metadat$Round[metadat$Bird.ID ==
areas_by_week$ID[x]][1])
areas_by_week$round <- paste("Round", areas_by_week$round)
call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
areas_by_week$treatment <- factor(areas_by_week$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
saveRDS(areas_by_week, "./data/processed/acoustic_space_area_by_individual_and_week.RDS")With rarefaction
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))
With rarefaction and compared to week 1
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {
X <- areas_by_week[areas_by_week$ID == x, ]
X$raref.area <- X$raref.area - X$raref.area[X$week == min(as.numeric(as.character(X$week)))]
return(X)
}))
ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))
Without rarefaction
ggplot(data = areas_by_week[areas_by_week$week != 5, ], aes(x = week, y = size, group = ID,
col = treatment)) + geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))
Without rarefaction and compared to week 1
areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {
X <- areas_by_week[areas_by_week$ID == x, ]
X$size <- X$size - X$size[X$week == min(as.numeric(as.character(X$week)))]
return(X)
}))
ggplot(data = areas_by_week, aes(x = week, y = size, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))
Acoustic space area by number of calls per individual and week
ggplot(data = areas_by_week, aes(x = n, y = size, col = treatment)) + geom_point() +
# geom_smooth(size = 1.2, method = 'lm') +
scale_color_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~round, ncol = 2, scales = "free_x") +
labs(y = "Acoustic space area", x = "Number of calls (n)") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.25))
- Acoustic space area by number of calls per individual and week by treatment
- Must have at least 2 calls
ggplot(data = areas_by_week, aes(x = n, y = log(size + 10), col = treatment)) + geom_point() +
geom_smooth(size = 1.2, method = "lm") + scale_color_viridis_d(begin = 0.1, end = 0.9) +
# facet_wrap(~ round, ncol = 2, scales = 'free_x') +
labs(y = "log(acoustic space area)", x = "Number of calls (n)") + theme_classic(base_size = 30) +
theme(legend.position = c(0.8, 0.7))
Acoustic space overlap
- 2 methods: pairwise distance between observations and kernel density overlap
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 2], ]
# run with density
dens_overlap_by_id <- space_similarity(X = Y, dimensions = c("PC1", "PC2"), group = "ID",
parallel = 10, pb = TRUE, outliers = 0.95, type = "mean.density.overlap")
dist_overlap_by_id <- space_similarity(X = Y, dimensions = c("PC1", "PC2"), group = "ID",
parallel = 10, pb = TRUE, outliers = 0.95, type = "distance")
saveRDS(list(dens_overlap_by_id = dens_overlap_by_id, dist_overlap_by_id = dist_overlap_by_id),
"./data/processed/acoustic_space_overlap_by_individual_2_metrics.RDS")Mantel test between the 2 methods
attach(readRDS("./data/processed/acoustic_space_overlap_by_individual_2_metrics.RDS"))
mantel.rtest(as.dist(1 - rectangular_to_triangular(dens_overlap_by_id, distance = FALSE)),
as.dist(rectangular_to_triangular(dist_overlap_by_id)), nrepet = 10000)## Monte-Carlo test
## Call: mantelnoneuclid(m1 = m1, m2 = m2, nrepet = nrepet)
##
## Observation: 0.6632477
##
## Based on 10000 replicates
## Simulated p-value: 9.999e-05
## Alternative hypothesis: greater
##
## Std.Obs Expectation Variance
## 5.5337964986 -0.0001518025 0.0143715608Individual overlap to its own group
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 2], ]
group_ovlp_l <- lapply(unique(Y$ID.week), function(x) {
# group data
X <- Y[Y$record.group == Y$record.group[Y$ID.week == x][1] & Y$week == Y$week[Y$ID.week ==
x][1], ]
# label all other individuals as group
X$ID2 <- ifelse(X$ID.week == x, "focal", "group")
# run with density
ovlp <- if (min(table(X$ID2)) > 4)
space_similarity(X = X, dimensions = c("PC1", "PC2"), group = "ID2", parallel = 10,
pb = FALSE, outliers = 0.95, type = "mean.mcp.overlap") else matrix(rep(NA, 3), nrow = 1)
dsts <- if (min(table(X$ID2)) > 1)
space_similarity(X = X, dimensions = c("PC1", "PC2"), group = "ID2", parallel = 10,
pb = FALSE, outliers = 0.95, type = "distance") else matrix(rep(NA, 3), nrow = 1)
out <- data.frame(ID = Y$ID[Y$ID.week == x][1], group = Y$record.group[Y$ID.week ==
x][1], week = Y$week[Y$ID.week == x][1], overlap.to.group = ovlp[1, 3], distance.to.group = dsts[1,
3], treatment = Y$treatment[Y$ID.week == x][1], round = Y$round[Y$ID.week ==
x][1])
return(out)
})
group_ovlp <- do.call(rbind, group_ovlp_l)
saveRDS(group_ovlp, "./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")
group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID), function(x) {
X <- group_ovlp[group_ovlp$ID == x, ]
X$overlap.to.group <- X$overlap.to.group - X$overlap.to.group[X$week == min(as.numeric(as.character(X$week)))]
X$distance.to.group <- X$distance.to.group - X$distance.to.group[X$week == min(as.numeric(as.character(X$week)))]
return(X)
}))
group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control", "Medium Stress",
"High Stress"))
ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to group acoustic space",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))
ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to group acoustic space",
x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
0.25))
# Partial mantel test between acoustic space
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID.week)
Y <- sels[sels$ID.week %in% names(tab)[tab >= 2], ]
treatment_labs <- sapply(unique(c(dens_overlap_by_id$group.1, dens_overlap_by_id$group.2)),
function(x) Y$treatment[Y$ID == x][1])
treatment_mat <- binary_matrix(labels = treatment_labs)
round_labs <- sapply(unique(c(dens_overlap_by_id$group.1, dens_overlap_by_id$group.2)),
function(x) Y$round[Y$ID == x][1])
round_mat <- binary_matrix(labels = round_labs)
# mantel.rtest(as.dist(treatment_mat), as.dist(1 - dist_overlap_by_id))
# mantel.partial(as.dist(treatment_mat), as.dist(1 - dens_overlap_by_id), zdis
# = as.dist(round_mat))sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID.week)
Y <- sels[sels$ID.week %in% names(tab)[tab >= 5], ]
# run with rarefaction
overlap_by_week <- space_similarity(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
parallel = 1, pb = TRUE, outliers = 0.95, type = "mean.mcp.overlap")
saveRDS(overlap_by_week, "./data/processed/acoustic_space_overlap_by_individual_and_week.RDS")trait_space_plot <- function(X, dimensions, group, indices, basecex = 1, title = "",
colors = c("#3E4A89FF", "#35B779FF"), point.alpha = 0.7, background.indices = NULL,
pch = 1, labels = c("sub-space", "total space"), legend = TRUE) {
total_coors <- as.ppp(as.matrix(X[, dimensions]), c(range(X[, dimensions[1]]),
range(X[, dimensions[2]])))
total_space <- raster(density.ppp(total_coors))
xlim <- range(X[, dimensions[1]])
ylim <- if (is.null(background.indices))
range(X[, dimensions[2]]) else c(min(X[, dimensions[2]]), c(max(X[, dimensions[2]]) + (max(X[, dimensions[2]]) -
min(X[, dimensions[2]])) * 0.11))
plot(x = X[, dimensions[1]], y = X[, dimensions[2]], col = "white", cex = basecex,
xlab = dimensions[1], ylab = dimensions[2], xlim = xlim, ylim = ylim, cex.lab = basecex,
xaxs = "i", yaxs = "i")
# add background group
if (!is.null(background.indices)) {
# get points
Y_bg <- X[background.indices, ]
if (nrow(Y_bg) > 1 & point.alpha > 0)
points(Y_bg[, dimensions[1]], Y_bg[, dimensions[2]], col = adjustcolor(colors[2],
alpha.f = point.alpha), pch = pch)
if (nrow(Y_bg) > 4) {
coors <- as.ppp(as.matrix(Y_bg[, dimensions]), c(range(Y_bg[, dimensions[1]]),
range(Y_bg[, dimensions[2]])))
raster_dens <- raster(density.ppp(coors))
palpre <- colorRampPalette(c("white", colors[2]))
colsn <- 10
palpre2 <- sapply(1:colsn, function(x) adjustcolor(col = palpre(colsn)[x],
alpha = seq(0.1, 0.9, length.out = 20)[x]))
image(raster_dens, add = TRUE, col = palpre2)
}
}
Y <- X[indices, ]
if (nrow(Y) > 1 & point.alpha > 0)
points(Y[, dimensions[1]], Y[, dimensions[2]], col = adjustcolor(colors[1],
alpha.f = point.alpha), pch = pch)
if (nrow(Y) > 4) {
coors <- as.ppp(as.matrix(Y[, dimensions]), c(range(Y[, dimensions[1]]),
range(Y[, dimensions[2]])))
raster_dens <- raster(density.ppp(coors))
palpre <- colorRampPalette(c("white", colors[1]))
colsn <- 10
palpre2 <- sapply(1:colsn, function(x) adjustcolor(col = palpre(colsn)[x],
alpha = seq(0.1, 0.9, length.out = 20)[x]))
image(raster_dens, add = TRUE, col = palpre2)
}
usr <- par()$usr
# legend
if (!is.null(background.indices) & legend)
legend(y = rep(usr[4] + ((usr[4] - usr[3]) * 0.05), 2), x = c(usr[1] + (usr[2] -
usr[1]) * 0.25, usr[1] + (usr[2] - usr[1]) * 0.75), col = colors, legend = labels,
pch = c(pch, pch), cex = basecex * 1.5, bty = "n", horiz = TRUE, xjust = 0)
title(title, cex.main = basecex * 1.2)
}
# trait_space_plot(X = sels, indices = which(sels$ID == '395WBHM'), dimensions
# = c('TSNE1', 'TSNE2'), point.alpha = 0.2, pch = 1) trait_space_plot(X = sels,
# indices = which(sels$ID == '395WBHM'), dimensions = c('TSNE1', 'TSNE2'),
# background.indices = which(sels$ID != '395WBHM' & sels$record.group ==
# sels$record.group[sels$ID == '395WBHM'][1]), point.alpha = 0.2, pch = 1)# Y <- sels[sels$ID == '395WBHM',] Y <- sels[sels$ID %in% c('395WBHM',
# '394WBHM', '360YGHM', '385YBHM'), ]
# trait_overlap(Y, dimensions = c('TSNE1', 'TSNE2'), group = 'ID', type =
# 'density', overlap.type = 'assymetric')
# X = Y dimensions = c('TSNE1', 'TSNE2') level <- '395WBHM' basecex <- 1 group
# <- 'ID' indices <- which(X[, group] == level)
trait_space_plot(X = sels, indices = which(sels$ID == "395WBHM"), dimensions = c("TSNE1",
"TSNE2"), background.indices = which(sels$ID != "395WBHM" & sels$record.group ==
sels$record.group[sels$ID == "395WBHM"][1]), point.alpha = 0.2, pch = 1)
# with points
out <- pblapply(unique(sels$ID), function(x) {
tiff(filename = file.path("./images", paste(x, sels$treatment[sels$ID == x][1],
"v2.tiff", sep = "-")), res = 120, width = 800, height = 800)
par(mfrow = c(2, 2))
# W <- sels[sels$ID == x, ]
for (i in 1:4) try(trait_space_plot(X = sels, dimensions = c("TSNE1", "TSNE2"),
indices = which(sels$ID == x & sels$week == i), background.indices = which(sels$ID !=
x & sels$week == i & sels$record.group == sels$record.group[sels$ID ==
x][1]), title = paste(x, "week", i), basecex = 1, labels = c(x, "group"),
legend = FALSE), silent = TRUE)
dev.off()
})Session information
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=pt_BR.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] PhenotypeSpace_0.1.0 warbleR_1.1.26 NatureSounds_1.0.4
## [4] seewave_2.1.8 tuneR_1.3.3.1 vegan_2.5-7
## [7] lattice_0.20-44 permute_0.9-7 raster_3.5-15
## [10] spatstat_2.3-3 spatstat.linnet_2.3-2 spatstat.core_2.4-0
## [13] rpart_4.1-15 nlme_3.1-152 spatstat.random_2.1-0
## [16] spatstat.geom_2.3-2 spatstat.data_2.1-2 GGally_2.1.2
## [19] adehabitatHR_0.4.19 adehabitatLT_0.3.25 CircStats_0.2-6
## [22] boot_1.3-28 adehabitatMA_0.3.14 ade4_1.7-18
## [25] deldir_1.0-6 sp_1.4-6 MASS_7.3-54
## [28] Rtsne_0.15 randomForest_4.7-1 formatR_1.11
## [31] knitr_1.37 kableExtra_1.3.4 ggplot2_3.3.5
## [34] viridis_0.6.2 viridisLite_0.4.0 pbapply_1.5-0
## [37] readxl_1.3.1 lubridate_1.8.0 remotes_2.4.2
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-3 rjson_0.2.21 ellipsis_0.3.2
## [4] rstudioapi_0.13 proxy_0.4-26 farver_2.1.0
## [7] fansi_1.0.2 xml2_1.3.3 codetools_0.2-18
## [10] splines_4.1.0 polyclip_1.10-0 jsonlite_1.8.0
## [13] cluster_2.1.2 spatstat.sparse_2.1-0 compiler_4.1.0
## [16] httr_1.4.2 assertthat_0.2.1 Matrix_1.3-4
## [19] fastmap_1.1.0 cli_3.2.0 htmltools_0.5.2
## [22] tools_4.1.0 gtable_0.3.0 glue_1.6.1
## [25] dplyr_1.0.8 Rcpp_1.0.8 cellranger_1.1.0
## [28] jquerylib_0.1.4 vctrs_0.3.8 svglite_2.1.0
## [31] xfun_0.29 stringr_1.4.0 rvest_1.0.2
## [34] lifecycle_1.0.1 goftest_1.2-3 terra_1.5-21
## [37] scales_1.1.1 spatstat.utils_2.3-0 parallel_4.1.0
## [40] RColorBrewer_1.1-2 yaml_2.3.5 gridExtra_2.3
## [43] sass_0.4.0 reshape_0.8.8 stringi_1.7.6
## [46] highr_0.9 rlang_1.0.1 pkgconfig_2.0.3
## [49] systemfonts_1.0.4 dtw_1.22-3 bitops_1.0-7
## [52] evaluate_0.15 purrr_0.3.4 tensor_1.5
## [55] labeling_0.4.2 tidyselect_1.1.2 plyr_1.8.6
## [58] magrittr_2.0.2 R6_2.5.1 fftw_1.0-6.1
## [61] generics_0.1.2 DBI_1.1.1 pillar_1.7.0
## [64] withr_2.4.3 mgcv_1.8-36 abind_1.4-5
## [67] RCurl_1.98-1.6 tibble_3.1.6 crayon_1.5.0
## [70] utf8_1.2.2 rmarkdown_2.11 grid_4.1.0
## [73] digest_0.6.29 webshot_0.5.2 signal_0.7-7
## [76] munsell_0.5.0 bslib_0.2.5.1